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The isoelectronic series of Be, Ne and Si are investigated using a variational determination of the 
second-order density matrix. A semidefinite program was developed that exploits all rotational and 
spin symmetries in the atomic system. We find that the method is capable of describing the strong 
static electron correlations due to the incipient degeneracy in the hydrogenic spectrum for increasing 
| central charge. Apart from the ground-state energy various other properties are extracted from the 

variationally determined second-order density matrix. The ionization energy is constructed using the 
extended Koopmans' theorem. The natural occupations are also studied, as well as the correlated 
hartree-Fock-like single particle energies. The exploitation of symmetry allows to study the basis 
set dependence and results are presented for correlation-consistent polarized valence double, triple 
and quadruple zeta basis sets. 
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I. INTRODUCTION 



The idea of a variational determination of the ground-state energy for a nonrelativistic many-body problem based 
, on the second-order density matrix (2DM) has a long history [J, 0, 0] and several highly appealing features. The 
energy of a system is a known linear functional of the 2DM. iV-particle wave functions never need to be manipulated 
since the energy is minimized directly in terms of the 2DM. However, the minimization is constrained because the 
variational search should be done exclusively with 2DMs that can be derived from an A^-particle wave function (or an 
ensemble of A^-particle wave functions). Such a 2DM is called ./V-representable, and the complexity of the many-body 
£S) . problem is in fact shifted to the characterization of this set of iV-representable 2DMs. The complete (necessary and 
£> 1 sufficient) set of conditions for iV-representability of a 2DM is not known in a constructive form, but it is clear that 
the energy from a minimization constrained by a set of necessary A^-representability conditions is a strict lower bound 
to the exact energy. Therefore this approach is highly complementary to the usual variational procedure based on a 
wave-function ansatz, which produces upper bounds. In addition the method is in principle exact, in the sense that 
. ■ as increasingly accurate set of ./V-representability conditions are imposed in the minimization, the resulting energy 
OS ; converges to the exact one. 

These are fascinating ideas for any true-blooded many-body theorist, as it comes close to the "ultimate reduction" 
of an interacting many-particle problem to solving a sequence of two-particle problems. In practice, however, imple- 
menting the method turns out to be very difficult and it is only in the last decade that serious attempts have been 
undertaken to turn the idea into a practical calculational scheme. The massive efforts by Mazziotti et al. [1, [H, @ and 
Nakata et al. [7], [8| are particularly notable. The main difficulty is of a technical nature: stringent iV-representability 
conditions require the positive semidefiniteness of matrix functionals of the 2DM, which turns the variational prob- 
lem into a so-called semidefinite program (SDP). Even applying the simplest "two-index" conditions, a direct energy 
minimization using Newton-Raphson methods requires a matrix operation scaling as M 12 (where M is the number of 
single-particle states) in each Newton-Raphson step. This can be circumvented in various ways, so that only matrix 
operations scaling as M e are needed. While these are nominally M e methods, the number of iterations required to 
reach convergence is very high and seems to rise with system size; in practice present implementations are probably 
about 100-1000 times slower than comparable methods such as CCSD. Still, one has the feeling that there is potential 
to turn it into a genuine M 6 method, and it is of interest to investigate the properties of SDP applied to various 
systems. 
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Up to now, most applications covered electronic structure calculations in atoms and molecules. Attention has 
been given primarily to the resulting energy. In this paper, we focus on three issues: (i) the performance of SDP in 
multireference situations (strong static correlations), (ii) the quality of the variationally obtained 2DM, and (iii) the 
dependence of the results on M (the size of the basis set). We do this by investigating three well-known examples 
in electronic structure theory: the isoelectronic series of Be, Ne, and Si. It is well known that the correlation energy 
for N electrons in the field of a positive point charge Z has a Z-dependence that strongly depends on N. For 
an increasing central charge Z, the hartree-Fock spectrum tends to the hydrogenic one, which has an "accidental" 
degeneracy related to a special symmetry in the Coulomb Hamiltonian. For the four-electron series the incipient 
degeneracy of the 2p and 2s orbitals leads to a vanishing particle-hole gap, inducing strong correlation effects with a 
correlation energy proportional to Z. For the ten-electron series this does not happen because a major shell is closed, 
and the correlation energy becomes flat for increasing Z. The 14-electron series again shows degeneracy effects, and 
in addition is a spin triplet. 

In Sec. II we provide theoretical and calculational backgrounds on the SDP implementation that is used. (Note 
that the techniques developed in this section have been used previously to study molecular dissociation 0, In 
particular, we pay attention to the way spin and rotational symmetry are exploited, enabling the use of quite large 
(cc-pVQZ) basis sets. In Sec. Ill the SDP results for the isolectronic series of Be, Ne, and Si are discussed. A summary 
is provided in Sec. IV. Atomic units are used throughout the paper. 



II. THEORY 



A. iV-representability conditions 

We will use second quantized notation where a' a (a a ) creates (annihilates) an electron in a single-particle (sp) state 
a. The Hamiltonian can be written as 

»7 afijd 

where t ai is the matrix element of the one-body part of the Hamiltonian (kinetic energy plus external potential) and 
V a p- n s is the antisymmetrized matrix element of the Coulomb interaction. The problem of finding the ground state 
of a quantum mechanical many-body system can be reformulated in terms of the second-order density matrix 

r a0nS = (^lalalasayl^) . (2) 

In principle T is a complex Hcrmitian matrix, but for a Coulomb Hamiltonian it is sufficient to consider real-symmetric 
matrices, 

r a( 3; 7 <5 = r 7< 5 ;Q/ 3 . (3) 

In addition T obeys the fermionic relations for antisymmetry in the sp indices 

r Q /3 ;7 5 = — Fp a - n s = — r a /3 ; 5 7 = r^Q,;^ . (4) 

The density matrix T can be determined variationally through the minimization of the energy functional 

E(T) = Tr [VH^] = \ r^rr^Sn* > (5) 

a/3;7<5 

where the reduced two-particle (tp) Hamiltonian is defined as 



H 



(2) 

The problem with this method is that the complete set of conditions that the density matrix has to fulfill to be 
derivable from a physical wave function (the so-called A^-representability conditions) is not known in a constructive 
form [111 ]. Therefore one minimizes the energy functional under a limited set of iV-representability conditions. Three 
simple conditions, known as the P, Q and G conditions 0,111] > are known to give quite good results. The P condition 
expresses the fact that the 2DM has to be positive semidefinite. The physical interpretation of the Q condition is 
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that the two-hole matrix, Q, has to be positive semidefinite; using basis anticommutation relations Q can be written 
as a homogeneous linear mapping, from the tp matrix space onto itself: 



) a ^s = {y N \a aaf3 ala\\y N ) 

= r aj g ;7 5 + — (S ay Sps — S a sSp 7 ) Tr r — Spspay + 5 a sPf3~/ ~ <5q 7 P/3<5 + <5/3 7 /O a< 5 



Here the particle number constraint has been used 

Tr T 

as well as the definition of the sp density matrix: 

Pa~f — 



N(N - 1) 



N 



Q/3;7/3 



(6) 



(7) 



The G condition demands that the particle-hole (ph) matrix G, is positive semidefinite; again, G can be written as a 
homogeneous linear mapping, from the tp matrix space onto the ph matrix space: 

G af3nS = (V N \aia ala 7 \y N ) 

= SpsPaj — r Q< 5 ;7 /5 . (8) 

Recently there has been progress on improved iV-representability conditions using the positive semidefiniteness of 
higher-order density matrices, e.g. the three-positivity conditions known as the Tj and T2 conditions [5I [T3|. Also 
some attempts have been made to improve ./V-representability while remaining strictly in tp space, by considering 
Hamiltonian dependent positivity conditions [Hj], or sharp bounds on the P, Q and G operators [3]. However, in 
the present paper we restrict ourselves to the standard P, Q and G conditions. 



B. Inclusion of spin symmetry 

1. General case 



When the Hamiltonian of the system is invariant under rotations in spin space, the eigenstates can be characterized 
by their total spin £ and spin projection /1. Explicitly introducing the electron spin, a sp state is written as = 
|asa)}, where a is the spatial orbital index and s a = ±5 is the spin projection. Two sp states can couple to a pair 
with total spin S = or S = 1. The corresponding pair creation operator is 

S 



B 



t 

ab;SM 



[4 <8> a\ 



A I 



s a s b 



ls a ±s b \SM)al Sa al s 



(9) 
(10) 



and the density matrix T in spin-coupled tp space is defined as 

E TgSf M ' = (^IflU^I^) • (") 
The B^B operator in Eq. (jTTJ) can now be further coupled to an object with good total spin. First one has to 
introduce B, 

Bcd-SM = (— l) S+JU -B c rf;S —Mi (12) 
which is again a good spherical tensor operator. Equation can now be rewritten as, 

I Si 



N I 



B 



ab:S 



B, 



cd;S' 



N \ 
E/Jt) 



(13) 



St 



The density matrices on the right of Eq. (| 13[) are classified by St = 0, 1, 2 and provide an equivalent representation 
of the 2DM of the /xth member of the spin multiplet. Note that the 2DMs of different members are trivially related 
through the Wigner-Eckart theorem, 



N I 



B 



ab;S 



Beds' 



S T 



I* 



[St] 

in terms of reduced matrix elements. Here [S] — ^J2S + 1. 



JV v 



B 



ab;S ' 



B, 



cd;S' 



(14) 
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2. Singlet ground state 

If the ground state has E — (spin singlet), the number of matrices involved in the minimization procedure is 
significantly reduced. Obviously for a singlet ground state the operator in Eq. (|13p has to be scalar, i.e., only the 
St = part is nonzero, and Eq. (|13[) reduces to 

00 r SM;S'M' r c pS f-i c\ 



where 



rf 6;cd = (^ol^SM^SM^o) ■ (16) 



is independent of M. This shows that, for a singlet ground state, the density matrix in a coupled tp basis is diagonal 
in S and M, and independent of the spin projection M. Instead of having to work with the full density matrix, 
all matrix manipulations can be performed on only two diagonal blocks, the S = and S = 1 matrices, which are 
respectively symmetric and antisymmetric in the indices related to the spatial orbitals. 

We now reformulate the minization problem in the spin-coupled representation. The Q-matrix in the coupled 
representation is similarly defined as: 

Q S a b ;c d = ^00\Ba b ;SMBi d .^M) • (17) 

It is clear that the Q matrix has an identical block diagonal structure as T. After some recoupling one can write the 
Q-mapping, from coupled tp space onto coupled tp space, as 

Qab;cd = r ffc;cd + -{Sachd + (-l) S S ad 5 bc )Tl T - p ac S bd - (- 1) S PbJad - Pbd^ac ~ (-1) S PadSbc , 



where the sp matrix p 



....<*M|aLac..l*m> (18) 
1 1 



2N-1^ 



■^>] 2 E r ^ as) 



and the trace 



Trr^E^E 1 ^* (20) 

S ab 



can be expressed in terms of the coupled r s . 

The G-matrix is a bit more involved. The coupled ph creation operator reads 



a Lsm = [4®ab] M (21) 



(-\)-2-^{\ Sa l- Sb \SM)al Sa a bsb 



The G-matrix in coupled ph space can now be written as: 

G S ab ;cd = «l4 b;S M^M|0 • (22) 

Again, one can prove that this matrix has the same block structure as the T and Q matrices. After some angular 
momentum recoupling we get the expression for the G-map in the coupled representation: 

G S ab:cd = hdPac - E \- S '^ |i 1 yj r ftf:c&- ( 23 ) 
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3. Nonsinglet states 

For higher-spin multiplets the same block decomposition is possible, provided a spin-averaged ensemble is considered. 
The density matrix for such an ensemble is defined in spin-coupled representation as 

E r£f = 2sTlE^I<SM^;S'M|^) . (24) 
n 

Note that the minimal energy can be reached for such an ensemble, since all members of the multiplet are degenerate. 
Performing the same manipulation as leading to Eq. (fT3")) and using the Wigner-Eckart theorem as in Eq. (fT4")l one 
obtains 

/ i\S'-M ( -ns— u 



ab,cd 



E [s]2 E< S ~ M|Sr0) ^frf- (S ^ " [Bl;s ® ||*£> • (25) 



Since (— l) s 'VP] = — /j,|00) one can use orthogonality of the Clebsch-Gordan coefficients to work out the sum 

over jj, in Eq. (f2"5j) . The result 



T,ySS':M = $SS' ,j r N\ 







#£) , (26) 



implies that the ensemble 2DM is again block-diagonal in spin, and the same formulas can be used as for the singlet 
case. 

C. Inclusion of rotational symmetry 

In atomic systems, the rotational symmetry of the Hamiltonian further reduces the dimension of the blocks involved 
in the density matrix. In exactly the same way as for spin, one can show that the density matrix of an ensemble, 
when avaraged over the third component of angular momentum, is diagonal in the two-particle angular momentum 
L and its z-component Ml, and completely independent of the value of Ml- What is more, for atomic systems there 
is also the parity (tt = ±1) of the two particle states. In the end one gets a density matrix that is composed out of 
blocks with fixed values for L n S, enabling one to solve the variational problem in large basis sets. The sp basis for 
systems with rotational and spin symmetry is written as \am a s a ), where a is shorthand for the radial basis state n a l a . 
The tp density matrix in spin and angular momentum coupled representation is defined as 

r£? = (* N \Bl h . iL „ s B cd ,L«s\* N ) , (27) 

where 

1 LS 
\M L M S 

= E E< z » m « l » m »i £M i)<^4 S6 ^ Ms,)o «- , -' n -'- ^ , »' n »''- ■ (28) 

m a mh s a Sb 

In an analogous way as for spin coupling, the spin and angular momentum coupled Q-matrix is defined as 

QiLT = (* N \B ab ^sBi d .^ s \* N ) , (29) 
out of which the coupled Q-map can be derived 

SuKiA± - tacS h lAtl > (30) 

with the sp density matrix defined as 

"Y^y _ l E ^2 r n a lal b l b ;n c l a n b l b • ( 31 ) 



^ab-L^S 



[4 ® 4 



2 2L 4 



(L*S) n b l b 
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The G-matrix is defined as 

= (**| [4*5*]™ (Hvai]™) 1 |*"> , (32) 
where again a is a spherical tensor operator defined as 

a bmbSb = (-l) lb+mb+ ^ +Sb a b - mb - Sb . (33) 
The spin and angular momentum coupled G-map from tp space on ph space becomes 



D. Energy optimization with a semidefinite program 

1. Interior point method 

The variational problem for the 2DM can be formulated as a so-called semidefinite program fl5| . a constrained 
optimization program where it is demanded that certain matrices, which are functions of the variables being optimized, 
remain positive semidefinite. In our case, there is a convex subspace of the matrix space, which is called the feasible 
region, where T, Q(T) and G(T) are positive semidefinite. In T-space, the direction of energy decrease is given by 
(— i?( 2 )) [see Eq. ([5])]. If the energy is to be minimized, the objective is to go as far as possible in this direction, without 
leaving the feasible region. The optimized density matrix is on the edge of the feasible region. Computationally, this 
problem is solved with an interior point method by optimizing the following cost function 

$(I\<) = Tr TH {2) -tlndetP(r) + G , (35) 

with 

/r o o \ 

T(T) = Q(T) . (36) 

\o o g(T)J 

The constant G has no influence on the solution but is added here in order to take into account the possibility that the 
matrices have certain explicit zero eigenvalues connected with imposing spin constraints (see the discussion following 
Eq. (|50[k in this case G can be considered an infinite constant). Starting from a large value of t, (e.g. t = 1), the cost 
function is minimized, and the resulting density matrix is used as a seed vector for the next minimization program 
with a smaller value of t. This procedure continues until convergence is reached for t — * 0, when the density matrix 
is at the edge of the feasible region. 



2. Implementation 

In addition to the positive semidefinite constraints, there are a number of linear constraints which the density 
matrix has to fulfill (e.g. particle number). These conditions are imposed by direct substitution. Suppose there are 
a number of linear constraints of the form, 

Tr TK [i) = fc w . (37) 

The way to impose these conditions is to limit the variations to the subspace orthogonal to the ifW's. Suppose the 
set of symmetric tp matrices {/ z } is an orthogonal basis of this subspace, 

Tr f l p = S zj Tr pK {3) = . (38) 
The tp density matrix can be expanded in the basis, 



r = E r ^+ c , 

i 



(39) 
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where C is a constant matrix obeying the inhomogeneous conditions 

Tr CK^ = k (l) . 



(40) 



For the minimization of the cost function at a certain value of t, Newton's method is used. At a given point Tq in 
matrix space, the gradient of the cost function is 



9$ 



= Tr fH^ - t [Tr f {T^ 1 + Q (Q(ro)- 1 ) + Q (S(ro)- 1 )}] 



(41) 



Using the Hermiticity of the Q and Q mappings [e.g. Tr Q(T)A = Tr Q(A)r], the gradient in matrix form reads: 



V<i> = E W-f = P [ H(2) - 1 ( r o 1 + 2 (QFo)' 1 ) + AS (^(ro)- 1 



(42) 



where P is the operator that projects onto the space spanned by the f l 's and A is the antisymmetrizer that projects 
ph space on tp space. The Hessian at To can be written as 



a 2 $ 

dTidTj 



t{Tr [/V/^ 1 ] +Tr [Q^QiTo^QU^Q^o)- 1 ] + Tr ^(/^(rorW^ro)- 1 ] } 



In Newton's method the search direction e is found by solving the linear system 



E 



a 2 $ 



9$ 



"r,<vr,'' ar< 



(43) 



(44) 



This system is solved using the linear conjugate gradient method [161 ]. In this method, only one matrix- vector 
multiplication is needed per iteration. The special structure of the Hessian can be exploited to construct a fast 
matrix-vector multiplication. The action of the Hessian on a tp matrix e is 



^tt^- = i{Tr [/T^eTo 1 ] + Tr [Q(f)Q(r )- 1 Q(e)Q(r )- 1 ] + Tr [Q{f )g(T )- 1 g(e)g(T )- 1 ] } 



(45) 



which can be written in matrix form as, 



He = tP 



rv^ro 1 + Q (Q(r r 1 Q(e)Q(T ) 



Ag (g(T )- 1 g(e)g(T Q y 



(46) 



It is clear that each conjugate gradient step can be calculated using only manipulations in the tp and ph matrix space. 

After the convergence of the conjugate gradient cycle, the direction of the Newton-Raphson step e is known. A line 
search in this direction is then performed in order to obtain the minimum of the cost function. Note that one always 
stays in the feasible region since the cost function goes to +oo at the edge. 



3. Imposing the spin constraints for E = 
The spin coupled form of the S z operator can be written as: 

* a 

This operator lives in ph-space, and we can force the vector 

{± z } S ab = ^SlSab , (48) 

to be an eigenvector of Q(T) with eigenvalue zero. In doing this we automatically impose the same constraints on 
G(T) for S x and S y due to the threefold degeneracy of the 5 = 1 block of the 2DM. It can be easily seen that in this 
case the expectation value of the total spin is zero, 

(^lE 2 !^) = (^lE 2 + ± 2 + £ 2 Z \^ N ) = . (49) 
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So the condition to be imposed on the density matrix becomes: 



s 



[S? 



1 1 



2N- 1 



(-If ? 2 



(50) 



For the projection of a tp density matrix on a spin singlet state, there are as many constraint matrices as there are sp 
matrix dimensions. Because of the zero eigenvalues in the G-matrix the projected density matrix is on the edge of the 
feasible region during the whole of the minimization process, and as a result, the cost function is infinity. This can be 
circumvented by taking the pseudo-inverse of the G-matrix, which excludes the E z -state from the inversion process. 
This will not alter the result of the program because the contribution of this state to the cost function is constant. 



4- Imposing the spin constraints for E / 

For higher-spin multiplets we use the spin-averaged ensemble (see Sec. Ill B 3|) . in which the 2DM has the same 
simple structure as for the singlet case. The expectation value of the E 2 spin operator is forced to be exact, using the 
linear constraint 

Tr r{£ 2 } = E(E + 1) , (51) 
where {E 2 } is the tp matrix representation of the E 2 operator, 



}ab:cd — 



(SacSbd + i-lfSadtbc) ■ (52) 



There is only one linear constraint for nonzero spin, in contrast to the numerous constraints for the projection onto a 
singlet state. It can therefore be expected that the spin constraints (i.e. the constraints on the 2DM ensuring that it 
is derivable from a wave function with good total spin) are less accurate than those for the singlet case. It is, in fact, 
known how to cure this situation [17] by considering not the spin-averaged ensemble but rather the 2DM derived from 
the highest-weight member (fi = E) of the multiplet. Similar to the spin singlet projection, one can then impose the 
condition that, since the spin-raising ladder operator E + destroys the wave function, the Q(T) matrix must have a 
zero eigenvalue (with an eigenvector in ph space corresponding to the £ + operator). In such a highest-weight schemes, 
the spin restrictions for the E ^ case are put on the same footing as for the singlet case; in fact, the highest-weight 
and the spin-averaged ensemble scheme are equivalent for the singlet case. However, the highest-weight scheme for 
E ^ requires one to keep track of more matrices and is computationally more demanding by about a factor of 10. 
We therefore used the ensemble scheme even for the nonsinglets (i.e. the Si atom), though we checked some cases with 
the highest-weight method for the spin. 



5. Spin and angular momentum projection 



When angular momentum is taken into account, everything becomes a bit more complicated, but the principles 
are the same as in the last paragraph. It can be shown that in a spin-and-angular-momentum-coupled basis the 
z-projections of E and A become 



1 V^rn r t - 
E z = -7=ZJI\ [ a nl® a m\ 

nl 

/2„ ,r , 1(1 + 0) 



(53) 
(54) 



Following the same argument as before, it can be imposed that the density matrix is derivable from an eigenstate 
with zero eigenvalue of respectively the E and A operators when 



(55) 
(56) 
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This can be translated into linear constraints on the 2DM, which are given in the Appendix. The projection on spin 
and angular momentum not equal to zero is again a less strict condition. The expectation values of A and £ are 
projected on the desired values: 



Tr T{£ 2 } = £(E+ 1) , 
Tr T{A 2 } = A(A + 1) , 

where the {S 2 } and {A 2 } are the tp matrix representations of the £ 2 and A 2 operators respectively. 



(57) 
(58) 



cd 



32- A 



2 A- 1 

x (S ac S bd + (-l) L+s+l " +l >S ad S bc ) , 



cd 



2- A 
A - 1 

x (S ac S bd + (-l) L+S+la+h S ad 5 bc ) 



(59) 



(60) 



III. RESULTS AND DISCUSSION 



Using the method explained in the previous section, the isoelectronic series of Be, Ne and Si were calculated from 
the neutral atom up to a central charge Z = 28. Beryllium and neon are both elements with a singlet ground state. In 
the silicon ground state the total spin and angular momentum are both one, which allows us to assess the quality of 
the spin and angular momentum constraints for S, A ^ 0. In order to study the basis set dependence, the properties 
of the ground state of the Be and Ne series were calculated in a cc-pVDZ a cc-pVTZ and a cc-pVQZ basis set [Tsl | . 
The Si series was only calculated in a cc-pVDZ and a cc-pVTZ basis set We used spherical harmonic (and not 
Cartesian) basis functions throughout. With the density matrices obtained from the SDP, several prop erties were 
studied. These are compared to estimates for non-relativistic energies based on experimental data 0, l2l| . and to the 
results of coupled cluster (CCSD) calculations, and in some cases, with full-configuration-interaction (CI) calculations. 

The basis functions used were those of the neutral atom, but with a rescaling r — ► rZ/N for the positive ions with 
Z > A. The CCSD and full-CI results were obtained using the MOLPRO program [Hj]. 



A. Ground-state energy 

The ground-state energies, calculated with various basis sets and methods, are shown in Tables HI HT1 and HTT1 for the 
Be, Ne, and Si isoelectronic series, respectively. Even in the best case (Be in cc-pVQZ), the calculated energies are at 
least 20 mhartree removed from the experimental estimate in [2(| Hl| . This is due to the difficulty of describing the 
interelectronic cusp in the exact wave function using finite sp basis sets. More relevant is the difference between 
the SDP (and CCSD) energies as compared to full-CI in the same basis set. This is shown in Figure Q] for the case 
of the Be series. Note that the CCSD energy is always above, the SDP energy below, the full-CI energy. For SDP, 
this simply reflects the nature of the variational problem. For the smallest cc-pVDZ basis set, SDP and CCSD have 
about the same level of accuracy. The difference with full-CI grows as the basis set size increases for both CCSD and 
SDP, but this effect is worse for the SDP. 

As far as the Z-dependence is concerned, the trend differs markedly for the cc-pV(D,T)Z and for the cc-pVQZ 
basis set. As Z increases there is a growing accuracy for the smaller basis sets in both CCSD and SDP, whereas for 
cc-pVQZ the accuracy decreases for CCSD and becomes constant for SDP. The reason for this difference is not clear, 
though it is probably connected to the incipient degeneracy of the 2s and 2p states and the quality of its description 
in the various basis sets, as is more fully described in the next Section. It should be noted that the SDP results are 
overall very accurate, even in the worst case (Z = 28, cc-pVQZ) differing less than 3 mhartree from full-CI. 

For the Ne series, full-CI calculations were only possible in the cc-pVDZ basis. From the results collected in Table [TT1 
it is seen that the SDP accuracy is significantly less than for Be, the largest deviation (28 mhartree) appearing for the 
neutral atom. This is actually consistent with PQG-condition SDP results for molecules, so it is likely that because 
of the small number of electrons the Be results are not representative. This is also borne out by the Si results in 
Table UTTl showing a maximal deviation between CCSD and SDP energies of 21 mhartree for the neutral atom. 



TABLE I: The ground-state energies of the Be series in the cc-pV(DTQ)Z basis sets using different methods. 





SDP 


HF 


4 


-14.617473 


-14.572338 


5 


-24.275712 


-24.216056 


6 


-36.387458 


-36.316267 


7 


-50.940925 


-50.860695 


8 


-67.931909 


-67.844323 


9 


-87.358767 


-87.265015 


10 


-109.22078 


-109.12175 


11 


-133.51761 


-133.414 


12 


-160.24908 


-160.14145 


13 


-189.41511 


-189.30392 


14 


-221.01564 


-220.90129 


15 


-255.05067 


-254.93347 


1 p. 

ID 


OQi s 


-zyl.4UU4 


17 


-330.42417 


-330.30206 


18 


-371.76264 


-371.63841 


19 


-415.5356 


-415.40942 


20 


-461.74304 


-461.61508 


21 


-510.38498 


-510.25537 


22 


-561.46143 


-561.3303 


23 


-614.97237 


-614.83983 


24 


-670.91783 


-670.78398 


25 


-729.29781 


-729.16274 


26 


-790.1123 


-789.97609 


27 


-853.36132 


-853.22404 


28 


-919.04486 


-918.90659 



cc-pVDZ 

CCSD 



full-CI 



SDP 



cc-pVTZ 
HF CCSD 



full-CI 



SDP 



cc-pVQZ basis 
HF CCSD 



expt. 



full-CI 



-24.27566 
-36.387421 
-50.940896 
-67.931884 
-87.358746 
-109.22076 
-133.51759 
-160.24906 
-189.41509 
-221.01563 
-255.05066 
-291.52017 
-330.42416 
-371.76263 
-415.53559 
-461.74303 
-510.38498 
-561.46142 
-614.97237 
-670.91783 
-729.2978 
-790.11229 
-853.36131 
-919.04486 



-14.61741 

-24.275684 

-36.387439 

-50.940909 

-67.931896 

-87.358755 

-109.22077 

-133.5176 

-160.24907 

-189.4151 

-221.01564 

-255.05066 

-291.52017 

-330.42416 

-371.76263 

-415.53559 

-461.74304 

-510.38498 

-561.46142 

-614.97237 

-670.91783 

-729.2978 

-790.1123 

-853.36131 

-919.04486 



-14.625431 
-24.300695 
-36.473162 
-51.137349 
-68.290965 
-87.932793 
-110.06209 
-134.67837 
-161.7813 
-191.37066 
-223.44627 
-258.00801 
-295.05582 
-334.58961 
-376.60935 

-421.115 
-468.10654 
-517.58394 
-569.5472 
-623.99631 
-680.93126 
-740.35204 
-802.25866 
-866.65111 
-933.5294 



-14.572873 
-24.234557 
-36.394215 
-51.045734 
-68.186797 
-87.816285 
-109.93353 
-134.53811 
-161.62969 
-191.20808 
-223.27309 
-257.82461 
-294.86255 
-334.38682 
-376.39737 
-420.89414 
-467.87712 
-517.34625 
-569.30152 
-623.7429 
-680.67038 



-14.623559 
-24.299207 
-36.471944 
-51.136311 
-68.290046 
-87.931958 
-110.06132 
-134.67764 
-161.78061 

-191.37 
-223.44563 
-258.0074 
-295.05522 
-334.58904 
-376.60879 
-421.11445 

-468.106 
-517.58342 
-569.54669 
-623.99581 
-680.93077 



-740.08394 -740.35156 
-801.98356 -802.25818 



-866.36925 
-933.24098 



-866.65064 
-933.52894 



-14.62381 
-24.29943 
-36.47214 
-51.136486 
-68.290206 
-87.932107 
-110.06146 
-134.67778 
-161.78074 
-191.37011 
-223.44574 
-258.00751 
-295.05533 
-334.58913 
-376.60888 
-421.11454 
-468.10609 
-517.5835 
-569.54677 
-623.99589 
-680.93084 
-740.35163 
-802.25825 
-866.65071 
-933.529 



-14.642807 
-24.321254 
-36.500934 
-51.177145 
-68.347448 
-88.010503 
-110.16555 
-134.81213 

-161.95 
-191.57899 

-223.699 
-258.30998 
-295.4119 
-335.00476 
-377.08857 
-421.66334 
-468.72912 
-518.28593 
-570.33383 
-624.87286 
-681.90309 
-741.42459 
-803.43742 
-867.94167 
-934.93744 



-14.572968 
-24.236385 
-36.40257 
-51.065945 
-68.22364 
-87.87405 
-110.01625 
-134.64967 
-161.77398 
-191.38895 
-223.49443 
-258.09031 
-295.17653 
-334.75303 
-376.81977 
-421.37673 
-468.42388 
-517.96121 
-569.9887 
-624.50634 
-681.51414 
-741.01206 
-803.00013 
-867.47832 
-934.44663 



-14.639589 
-24.317643 
-36.497178 
-51.173335 
-68.343621 
-88.00667 
-110.16171 
-134.8083 
-161.94616 
-191.57514 
-223.69514 
-258.3061 
-295.40801 
-335.00085 
-377.08463 
-421.65938 
-468.72513 
-518.28191 
-570.32977 
-624.86877 
-681.89895 
-741.4204 
-803.43318 
-867.93738 
-934.93309 



-14.640124 
-24.31822 
-36.497761 
-51.173918 
-68.344203 
-88.007254 
-110.1623 
-134.80889 
-161.94677 
-191.57575 
-223.69577 
-258.30675 
-295.40867 
-335.00153 
-377.08533 
-421.66011 
-468.72588 
-518.28269 
-570.33058 
-624.86961 
-681.89983 
-741.42132 
-803.43414 
-867.93838 
-934.93413 



-14.66736 
-24.34892 
-36.53493 
-51.22284 
-68.41171 
-88.10113 
-110.29089 
-134.98088 
-162.17102 
-191.86127 
-224.0516 
-258.742 
-295.93244 
-335.62293 
-377.81344 
-422.50398 
-469.69455 
-519.38513 
-571.57572 
-626.26633 
-683.45695 
-743.14758 
-805.33822 
-870.02886 
-937.21951 



TABLE II: The ground-state energies of the Ne series in the cc-pV(DTQ)Z basis sets using different methods. 



z 


SDP 


cc- 

HF 


pVDZ 

CCSD 


full-CI 


SDP 


cc-pVTZ 
HF 


CCSD 


SDP 


cc-pVQZ basis 
HF 


CCSD 


expt. 


10 


-128.70843 


-128.48878 


-128.67964 


-128.68088 


-128.86088 


-128.53186 


-128.81081 


-128.92686 


-128.54347 


-128.87106 


-128.9376 


11 


-161.80049 


-161.59591 


-161.77283 


-161.77411 


-161.97703 


-161.65496 


-161.92829 


-162.05038 


-161.67155 


-161.99595 


-162.0659 


12 


-198.88784 


-198.70208 


-198.86199 


-198.86309 


-199.11372 


-198.79861 


-199.06598 


-199.19913 


-198.82303 


-199.14502 


-199.2204 


13 


-239.97194 


-239.80393 


-239.94802 


-239.94883 


-240.26728 


-239.9582 


-240.22028 


-240.36392 


-239.98965 


-240.30977 


-240.3914 


14 


-285.04223 


-284.88894 


-285.02004 


-285.02061 


-285.43166 


-285.12786 


-285.38525 


-285.53886 


-285.16605 


-285.48453 


-285.5738 


15 


-334.08381 


-333.94195 


-334.06299 


-334.06338 


-334.6021 


-334.30313 


-334.55624 


-334.72067 


-334.3492 


-334.66615 


-334.7642 


16 


-387.08194 


-386.94882 


-387.06219 


-387.06246 


-387.77553 


-387.48107 


-387.73017 


-387.90757 


-387.53731 


-387.85281 


-387.9608 


17 


-444.02427 


-443.89781 


-444.00531 


-444.00551 


-444.95002 


-444.6598 


-444.9051 


-445.09826 


-444.72921 


-445.04331 


-445.1622 


18 


-504.90101 


-504.77972 


-504.88268 


-504.88282 


-506.12426 


-505.83808 


-506.07977 


-506.29194 


-505.92402 


-506.23681 


-506.3673 


19 


-569.70474 


-569.58754 


-569.68689 


-569.68701 


-571.29734 


-571.01502 


-571.25329 


-571.48788 


-571.12105 


-571.43261 


-571.5754 


20 


-638.42981 


-638.31595 


-638.41239 


-638.41248 


-640.46864 


-640.18995 


-640.42496 


-640.68561 


-640.31973 


-640.63013 


-640.7891 


21 


-711.07205 


-710.96091 


-711.05497 


-711.05505 


-713.63756 


-713.36231 


-713.59424 


-713.88444 


-713.51958 


-713.8289 


-713.9988 


22 


-787.6282 


-787.51937 


-787.61143 


-787.6115 


-790.80365 


-790.53161 


-790.76063 


-791.08414 


-790.72018 


-791.02848 


-791.2132 


23 


-868.09589 


-867.98895 


-868.07932 


-868.07938 


-871.96637 


-871.69743 


-871.9237 


-872.28418 


-871.92117 


-872.22853 


-872.4291 


24 


-952.47304 


-952.36781 


-952.45674 


-952.45679 


-957.12544 


-956.85936 


-957.08305 


-957.48456 


-957.12224 


-957.42872 


-957.6463 


25 


-1040.7584 


-1040.6545 


-1040.7422 


-1040.7422 


-1046.2804 


-1046.0171 


-1046.2383 


-1046.6846 


-1046.3231 


-1046.6288 


-1046.8646 


26 


-1132.9505 


-1132.8479 


-1132.9345 


-1132.9345 


-1139.431 


-1139.1702 


-1139.3892 


-1139.8845 


-1139.5236 


-1139.8285 


-1140.0838 


27 


-1229.0485 


-1228.9471 


-1229.0327 


-1229.0328 


-1236.5769 


-1236.3184 


-1236.5353 


-1237.0837 


-1236.7235 


-1237.0277 


-1237.3039 


28 


-1329.0518 


-1328.9514 


-1329.0361 


-1329.0361 


-1337.7178 


-1337.4616 


-1337.6764 


-1338.2821 


-1337.9226 


-1338.2261 


-1338.5247 



TABLE III: The ground-state energies of the Si series in the cc-pV(DT)Z basis sets using different methods. The results under SDP were calculated using the ensemble 
avaraged spin projection, those under SDP* were calculated using the maximal weight method. 



z 


SDP 


SDP* 


cc-pVDZ 

HF 


CCSD 


SDP 


cc-pVTZ 
HF 


CCSD 


expt. 


14 


-288.93962 


-288.92921 


-288.84644 


-288.91895 


-289.02515 


-288.85215 


-288.9835 


-289.359 


15 


-340.36765 




-340.27338 


-340.34709 


-340.50472 


-340.33467 


-340.46205 


-340.872 


16 


-396.10801 




-396.01679 


-396.08749 


-396.43974 


-396.27384 


-396.39711 


-396.869 


17 


-456.09635 




-456.00926 


-456.0759 


-456.80372 


-456.64236 


-456.7617 


-457.337 


18 


-520.29362 


-520.27860 


-520.21067 


-520.27348 


-521.58 


-521.42294 


-521.53858 


-522.269 


19 


-588.68067 




-588.60149 


-588.66097 


-590.75683 


-590.60373 


-590.7159 


-591.66 


20 


-661.24791 




-661.17202 


-661.22871 


-664.32396 


-664.17433 


-664.28328 


-665.507 


21 


-737.99017 




-737.91714 


-737.97148 


-742.2714 


-742.12467 


-742.23072 


-743.808 


22 


-818.90446 


-818.88808 


-818.83388 


-818.88621 


-824.58949 


-824.44532 


-824.54882 


-826.559 


23 


-903.98889 




-903.92042 


-903.97105 


-911.27007 


-911.12778 


-911.22912 


-913.762 


24 


-993.24222 




-993.1756 


-993.22475 


-1002.3054 


-1002.1648 


-1002.2643 


-1005.413 


25 


-1086.6636 




-1086.5986 


-1086.6465 


-1097.6895 


-1097.5501 


-1097.6482 


-1101.513 


26 


-1184.2525 


-1184.2381 


-1184.1889 


-1184.2357 


-1197.4173 


-1197.2789 


-1197.3759 


-1202.061 


27 


-1286.0084 




-1285.9461 


-1285.9919 


-1301.4847 


-1301.347 


-1301.4431 


-1307.057 


28 


-1391.9311 




-1391.8699 


-1391.9147 


-1409.8882 


-1409.7512 


-1409.8466 


-1416.5 
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FIG. 1: Difference between approximate (CCSD or SDP) and full-CI energies for the Be series in all three basis sets. 



B. Correlation energy 



The correlation energies were calculated by taking the difference of the SDP and CCSD energies with the hartree- 
Fock results in the same basis set. The results labeled "experimental" are the estimates in [2l| . 



1. Beryllium series 



In Fig.[5]the SDP correlation energy is shown as a function of central charge Z for the different basis sets. Note that 
on the plot the difference between the CCSD and full-CI correlation energies would not be visible. The experimental 
curve is linear in Z, as a direct consequence of the near-degeneracy of the ground state [2l|. One can calculate a 
perturbative series expansion of the exact and hartree-Fock energy in powers of the corresponding series for the 
correlation energy starts with a constant if the hydrogenic ground state is nondegenerate, or with a linear term in Z in 
case of degeneracy. The SDP correlation energy does not follow this trend: it goes linear in the beginning, but becomes 
concave in the cc-pVDZ and cc-pVTZ basis, or convex in the cc-pVQZ basis. This failure, however, is not related to 
the SDP method as the trend is the same in full-CI. It simply reflects the fact that the incipient degeneracy is not well 
described in these basis sets. This can also be seen by calculating the Z = 1 hydrogen spectrum (corresponding to the 
Z — > oo situation, when the electron-electron interaction can be neglected) in the basis sets: the 2s and 2p energies 
are not degenerate, but differ by 5.8 mhartree (cc-pVDZ), 2.0 mhartree (cc-pVTZ) and -2.3 mhartree (cc-pVQZ). 
Note that for CC-pVQZ the 2p energy actually drops below the 2s energy, explaining the different (convex/concave) 
behavior of the curves. To make sure we also performed calculations in the cc-pVDZ basis after rescaling (r — > ar) 
it in such a way that the hydrogenic 2s-2p degeneracy is exact. In this basis the SDP correlation energy (also shown 
in Fig. [2]) indeed has the correct linear behavior. It is clear from the above discussion that SDP is indeed capable of 
providing accurate correlation energies in the presence of near degeneracies, when other many-body techniques (like 
density functional theory or MP2) can fail. 



2. Neon series 



In Fig. [3] the correlation energy is shown for all three basis sets as a function of Z. Because Ne is a closed shell 
atom, there is no near-degeneracy for large Z values and the exact correlation should be asymptotically constant in Z, 
as is indeed visible in the experimental curve. Due to basis set effects, this constant behavior is imperfectly realized, 
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FIG. 2: The SDP correlation energy for the Be series in all three basis sets, and in a rescaled basis set that exhibits hydrogen-like 
behaviour (degeneracy between the 2s and 2p level). For comparison, the CCSD and experimental values are also shown. 

but the SDP follows the same trends as CCSD for all basis sets. Note that the approximation to a constant behavior 
at large Z is best for the largest basis set. The decrease in correlation energy for increasing Z, in contrast to the 
slight rise in the experimental correlation energy, can be attributed to the fact that the basis sets were optimized for 
the neutral atom. While the rescaling procedure fixes the nuclear cusp, the resulting basis set is obviously far from 
optimal for highly charged ions. 

3. Silicon series 

For silicon, only the cc-pVDZ and cc-pVTZ basis have been used [Fig. @]. As was the case for Be, the theoretical 
linear rise with Z is thwarted by imperfections in the basis sets. However the SDP correlation energy closely tracks 
the CCSD one. The Si ground state is a spin triplet. The results in Table Hill have been obtained using the spin- 
averaged ensemble, as explained in Sec II. In order to assess the quality of the spin constraints, we have also performed 
calculations using the highest-weight method, for Z=14, 18, 22, and 26 with the cc-pVDZ basis set. The resulting 
energies are also reported in Table IIII1 The energy differences between the approaches are sizeable, with differences 
as large as 20 mhartree, reflecting the weaker nature of the spin constraints imposed in the spin-averaged scheme. 
However, the discrepancy between the two approaches is stable for increasing Z. 

C. Ionization energies 

Other properties can be used to gauge the quality of the 2DM, e.g., the ionization energies of the different atomic 
ions, which can be easily calculated using the extended Koopmans' theorem (EKT) [23|, 0, [25| . The EKT provides 
a single particle picture of the ground state, with sp energies and spectroscopic factors. The ionization energies are 
shown in Figure the agreement between calculated and experimental values is very good, pointing to the realistic 
nature of the variationally obtained 2DM. The good agreement with experiment reflects the fact that the error in the 
description of the interelectronic cusp largely cancels since the ionization energy is an energy difference. For Be and 
Ne it is clear that the basis set limit is nearly reached at the cc-pVTZ - cc-pVQZ level. Even for Si the experimental 
ionization energy is closely reproduced. 



15 



0.4 - 



0.35 




V V V 



V V V 



V V V v 



V V V V V 



o O o o 



0.3 - 



o O o o o 



o O O o o 



0.25 



1 0.2 



n □ 



□ □ 



□ □ 



□ □ 



0.15 



0.1 



0.05 



DZ SDP + 
DZ CCSD x 
TZ SDP A 
TZ CCSD □ 
QZ SDP v 
QZCCSD o 
expt. • 



t : + + + + + + + 



X X X x 



10 



15 



20 



25 



FIG. 3: The SDP correlation energy for the Ne series in all three basis sets. For comparison, the CCSD and experimental 
values are also shown. 



D. Correlated hartree-Fock-like single-particle energies 



A different sp picture is given by the correlated hartree-Fock-like sp orbitals and energies. These are constructed 
by diagonalizing the sp Hamiltonian: 

h ai = (T + U) ai + ^2 v ap- n sP(iS , (61) 

08 

where the first-order density matrix (1DM) (p) is constructed from the variationally determined 2DM. As an example 
of this method, the sp energies for the isoelectronic series of Be in a cc-pVDZ basis are shown in Figure [H Notice 
that when Z increases, the energy levels approach those of the hydrogen atom. Similar behavior is present for the 
other basis sets and for the Ne and Si isoelectronic series. 



E. Natural occupations 



The eigenvalues of the sp density matrix (i.e. the natural orbital occupation numbers) provide insight into the 
extent of correlation. The occupation numbers from SDP are always very close to those from full-CI, differing by at 
most 0.005. Of particular interest are the occupations of the quasidegenerate 2s and 2p orbitals in Be. These are 
shown in Fig. [7] The sum of the 2s and 2p occupations is nearly 1 and increasingly so for large Z. This implies that 
only the 2s and 2p are partially occupied in the large- Z limit. The shapes of the curves reflect the aforementioned 
imperfections in the basis sets, with the 2s below the 2p for cc-pVDZ and cc-pVTZ, and above the 2p for cc-pVQZ. 



IV. SUMMARY 



Variational methods based on the second-order density matrix seem to hold great promise as an ab initio many-body 
technique, but there is room for improvement , especially as regards computational efficiency (improved algorithms) 
and accuracy (better characterization of the Af-representable set). We investigated the isoelectronic series of Be, Ne 
and Si using the P, Q and G A^-representability conditions. A significant speedup is obtained when spin and rotational 
symmetry is taken into account. This allowed us to investigate the properties of the SDP method with increasing basis 
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set size (cc-pV(D,T,Q)Z). The energies so obtained are reasonably accurate, but the accuracy seems to diminish with 
increasing basis set size. The SDP method is capable of describing the strong static electron correlations appearing 
in the beryllium and silicon series due to the incipient degeneracy in the hydrogenic spectrum for increasing central 
charge. The ionization energies, constructed using the extended Koopmans' theorem, are surprisingly good. Also the 
natural occupations are reproduced very well when compared to full-CI results in the same basis sets. Hence, the 
physical content of the variationally determined second-order density matrix seems to be reliable. 

Apart from a study of the potential energy surface for some diatomic 14-electron molecules [Io|, we intend to 
investigate fermionic and bosonic Hubbard models on one and two - dimensional lattices. Further work is also needed 
to ameliorate the computational cost of the method, and to increase the accuracy without introducing three-index 
conditions. 



APPENDIX: CONSTRAINT MATRICES 



The linear constraints for imposing the spin singlet condition are given by: 

V k < I : Tr T [kl] K = , 
where the constraint matrices have the following form 

[kl] tsS [kl] f S ,( i\S[kl] f S j_l i\S[kl]fS ,[kl] f S 

K ab;cd = Jab;cd+{- 1 ) Jba;cd + Jab;dc + Jba;dc 



,[kl] f S , / I ^[kl] f S , /_■■ ^[kl] rS , [kl] fS 
+ Jcd;ab'\ L J Jdc-ab^\ 1 ) J cd-,ba ' Jdc;ba ' 



with 



[kl] fS 



J ab:cd 



1 1 



-(-1) 




2N - 1 

The constraints for spin and angular momentum singlet projection are: 

V k<l : Tr T [ j l] K = , 
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FIG. 5: The ionization energy scaled with for the Be, Ne and Si series in the different basis sets compared with experimental 
results. 
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FIG. 6: The single particle levels obtained in a correlated hartree-Fock-like scheme (see text) for the Be series in a cc-pVDZ 
basis set. 



where the constraint matrices have the following form 

[kl] AL^S) ( _,-,L+S+l a +l b [kl] fl 
JJab;cd \ 1 ) J-lb 

, [kl] rL'S ,( i\S+L+l a +l b [kl] fL 'S , [kl] rL'S , / i \L+S+l a +l b [kl] f L w S /a c\ 

"+" jJcd-ab T" I L ) jJdc;ab + jJdc;ba "+" I X J jJcd-ba > l A -°; 



[fei] ^(L"S) _ [kl] f (L*S) , _ U L+S+I a +I b [kl] rL'S , [kl] f (L"S) , _ u S+L+l a +l b [kl] f (L"S) 
J^ab-cd — JJab-cd ' \ 1 ) J J ba;cd ' J J ba;dc ' I 1 1 J -I ab;dc 
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FIG. 7: The natural occupation of the 2p-orbital and one minus the occupation of 2s-orbital for the Be series, in all three basis 
sets. 



where J can mean either E or A and 
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